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ABSTRACT 

O 

o 

CN ' We present 3D magnetohydrodynamic (MHD) numerical simulations of the 

evolution of self-gravitating and weakly magnetized disks with an adiabatic equa- 
C/^ ■ tion of state. Such disks are subject to the development of both the magnetoro- 

\Q • tational and gravitational instabilities, which transport angular momentum out- 

ward. As in previous studies, our hydrodynamical simulations show the growth 
of strong m = 2 spiral structure. This spiral disturbance drives matter toward 
^ ' the central object and disappears when the Toomre parameter Q has increased 

O ■ well above unity. When a weak magnetic field is present as well, the magnetoro- 

Q\ • tational instability grows and leads to turbulence. In that case, the strength of 

O 

O ■ hydrodynamical run and oscillates periodically, reaching very small values at its 

• minimum. We attribute this behavior to the presence of a second spiral mode 

with higher pattern speed than the one which dominates in the hydrodynamical 
2 ■ simulations. It is apparently excited by the high frequency motions associated 

^ ■ with MHD turbulence. The nonlinear coupling between these two spiral modes 

gives rise to a stress tensor that oscillates with a frequency which is a combination 



the gravitational stress tensor is lowered by a factor of about 2 compared to the 



of the frequencies of each of the modes. This interaction between MHD turbu- 
lence and gravitational instabilities therefore results in a smaller mass accretion 
rate onto the central object. 
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1. Introduction 



In systems such as the disks surrounding low mass protostars or active galactic nuclei, 
the simultaneous appearance of both gravitational and magnetic instabilities is expected. 
During the first stages of their evolution, for example, protoplanetary disks are expected 
to be rather massive because of strong infall from the parent molecular cloud. As the disk 
builds up in mass as a result of the collapse of an envelope, its surface mass density becomes 
large enough for gravitational instabihties to develop (e.g., Laughlin & Bodenheimer 1994). 
These disks are also believed to be sufficiently ionized, at least over some extended regions, 
to be coupled to a magnetic field (Gammie 1996; Sano et al. 2000; Fromang et al. 2002). 

By modeling the outer parts of disks around quasi-stellar objects (QSOs) as steady, 
viscous, geometrically thin, and optically thick, Goodman (2003) has argued that they are 
self-gravitating. More precisely, he predicts self-gravitational instabilities to develop beyond 
about 10~^ parsecs from the central object. In addition, it has been suggested by Menou & 
Quataert (2001) that self-gravitating regions of disks around QSOs are likely to be coupled 
to a magnetic field. 

The stability of a thin, self-gravitating gas disk is controlled by the Toomre Q parameter 
(Toomre 1964): 



where Cs is the sound speed, k is the epicyclic frequency (see, e.g., Binney & Tremaine 1987), 
S is the disk surface mass density and G is the gravitational constant. Gaseous disks are 
unstable against axisymmetric perturbations when Q < 1, and against non-axisymmetric 
perturbations when Q^l. 

Since analytical predictions of the nonlinear evolution of gravitational instabilities are 
difficult, there have been a large number of numerical simulations of gravitationally unsta- 
ble disks. Despite the rather daunting technical problems of combining three-dimensional 
(3D) hydrodynamic calculations with rapid and accurate Poisson equation solvers, signifi- 
cant progress have been made. To do so, the energetics must be treated crudely, with the 
focus squarely on purely dynamical behavior. Using this strategy, the above Q criterion for 
instability has been confirmed (and shown to still be approximately valid for disks of finite 
thickness), and the properties of the unstable modes have been studied as a function of the 
disk parameters (Tohline & Hachisu 1990; Woodward et al. 1994). Several authors have 
investigated the saturation properties of the instability, and have shown that it is capable 
of transporting significant amount of mass and angular momentum in a few orbital times 
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(Papaloizou & Savonije 1991; Laughlin & Bodenheimer 1994; Pickett et al. 1996; Laughlin 
et al. 1997). The first calculations mostly used simple adiabatic equations of state (EOS). 
More recently, isothermal disks have also been studied (Pickett et al. 1998, 2000; Boss 1998; 
Mayer et al. 2002). Some new investigations also include a simplified treatment of the disk 
radiative coohng (Pickett et al. 2003; Rice et al. 2003; Boss 2002). 

All these models were purely hydrodynamical, and neglected the effect of magnetic fields. 
However, it is known that stability of astrophysical disks is extremely sensitive to the presence 
of weak magnetic fields. In particular, the magnetorotational instability (MRI) completely 
disrupts laminar Kcplcrian fiow when a subthcrmal magnetic field of any geometry is present. 
This was first understood by Balbus & Hawley (1991). Since then, it has been shown through 
many numerical simulations that the nonlinear outcome of the MRI is MHD turbulence, 
which, in common with gravitational instabilities, transports angular momentum outward 
(see Balbus & Hawley 1998, or Balbus 2003, for a review). Since disks around low-mass 
stars and around QSOs may be both magnetized and self-gravitating, the spiral structure 
gravitational transport described above must somehow develop in a medium in the throes 
of MHD turbulence. 

The question naturally arises as to how these two powerful instabilities interact with 
one another. What is the ultimate effect on the global properties of accretion disks, and in 
particular, on the critical transport properties of mass and angular momentum? To keep 
this initial investigation tractable, we must restrict ourselves here to an adiabatic EOS. But 
the dynamical behavior of "simple" adiabatic disks is still rich, and contains unanticipated 
findings. In a companion paper to this one (Fromang et al. 2004, hereafter paper I), we carried 
out 2D axisymmetric numerical simulations of the evolution of massive and magnetized disks. 
The results show that the MRI behaves in a self-gravitating environment as it does in zero 
mass disks. Turbulent transport of angular momentum causes the disk to evolve toward 
a two component structure: (1) an inner thin disk in Keplerian rotation fed by (2) an 
outer thick disk whose rotation profile deviates from Keplerian, strongly infiuenced by self- 
gravity. However, angular momentum transport by gravitational instabilities cannot develop 
in axisymmetric simulations, which leaves unanswered the question of the outcome of the 
interaction between both instabilities. This is the subject of the present paper. 

The plan of the paper is as follows: in section 2, we present our numerical methods. 
The initial state of our simulations will be described in section 3. We present our results in 
section 4 and, finally, give our conclusions in section 5. 
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2. Numerical methods 
2.1. Algorithms 

The calculations in this paper are based on the equations of ideal MHD: 



^ + V.(pv) = 0, (2) 

/(9v \ 1 

pNt + v-Vv =-VP-pV$ + — (VxB)xB, (3) 

\ot J 47r 

— = Vx(vxB), (5) 



where p is the mass density, e is the energy density, v is the fluid velocity, B is the magnetic 
field, P is the gas pressure and $ = $s + $c is the total gravitational potential, which 
has contributions $s from the disk self-gravity and $c from a central mass. The Poisson 
equation determines the gravitational potential, 

V2$, = AttGp, (6) 

and to close our system of equations, we adopt an adiabatic equation of state for a monoatomic 
gas: 

P=(7-l)e, 7 = 5/3. (7) 

To solve these equations, we use the GLOBAL code (Hawley & Stone 1995). This uses 
standard cylindrical coordinates (r, 0, z) and time-explicit Eulerian finite differences. The 
magnetic field is evolved using the combined Method of Characteristics and Constrained 
Transport algorithm (MOC-CT), which preserves the divergence of the magnetic field to 
machine accuracy. Finally, we use outfiow boundary conditions in the radial and vertical 
directions, and periodic boundary conditions in (p. 

In its original form, GLOBAL did not include a Poisson solver, and the development of 
such a routine represents a major technical component of the results we report here. The 
calculation is done in two steps. The potential $5 is first computed at the grid boundary, 
using the spectral decomposition decribed below, and then calculated on the whole grid using 
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a very rapid method. It is the first step, the boundary calculation, that is computationally 
expensive. 

In the expansion of $5, we have adopted the method of Cohl & Tohline (1999), which 
uses half-integer Legendre functions in the Green's function. This method is better suited 
to cylindrical coordinates than the traditional expansion in spherical harmonics, which are 
of course tailored to spherical coordinates. Following Cohl & Tohline (1999), $s may be 
written 



$,(r, 0, ^) = — ^ / dr' P^^'^fj''^ emQm-i/2{x) cosm(</. - 0') . (8) 

Here, dr' = r'dr'd(f)'dz' is the elementary volume element, and the integral is taken over the 
whole computational domain. Qm-1/2 denotes the half-integer order Legendre function of 
the second type (Abramowitz & Stegun 1965). The argument x is a function of position: 

^^r2^r^^^ (9) 

The Legendre functions are computed once at the beginning of each simulation and stored 
in memory. At each time step, we calculate $s using equation (8), in which the sum over m 
is truncated at some upper value rrimax- We then calculate $s everywhere on the grid, using 
a combination of a Fourier transform in and the 2D Successive Over Relaxation (SOR) 
Method (Hirsch 1988) in the (r, z) plane. Although this is an efficient method, the calculation 
of the self-gravitating potential is still very demanding of computational resources. For the 
resolution (A^^j A^(/), A^z) = (128, 64, 128) used in this paper, the time required by the Poisson 
solver still represents ~ 40% of the computation time for rrimax = 8. 



2.2. Diagnostics 

We introduce and define some key quantities that have been used to analyze the results 
of the simulations. We denote the ratio of the volume averaged thermal pressure to the 
volume averaged magnetic pressure as (/?): 



This parameter is used primarily as a measure of the initial magnetic field strength. 
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In 3D numerical simulations of magnetized self-gravitating disks, angular momentum is 

transported by the sum of the Maxwell, Reynolds, and gravitational stress tensors. Following 
Balbus & Papaloizou (1999) and Hawley (2000), we define the height and azimuthal averages 
(noted with an overbar) of each these respective stresses as: 



TrT{r,t) = (11) 
T^f-{r,t) = pivW-^^, (12) 

As in paper I, volume averages of these quantities will be denoted as {T^"'^){t), etc. Note 
that Tj^™^ is associated with the gravitational torque resulting from non-axisymmetric disk 
structure. This quantity clear vanishes in an axisymmctric simulation {rrimax = 0)- In this 
case, the standard a parameter (Shakura & Sunyaev 1973) can be defined as the sum of the 
Maxwell and Reynolds stress tensors normalized by the gas pressure: 

a{r,t) ^ — . (14) 

P{r,t) 



3. Initial model 



We start our simulations with a disk model which is as close as possible to hydrostatic 
equilibrium: 



-VP - pV ($, + $c) + prn'^Br = . (15) 

Here Q is the angular velocity and e^. is the unit vector in the radial direction. The coordinate 
system has its origin on the disk center. The potential $c is due to a central mass Mp. We 
chose Mc — 2Md, where is the disk mass. The initial disk model is gravitationally 
unstable. 

Because of the presence of the disk self-gravity, equation (15) has to be solved iteratively. 
We use the Self-Consistent Field (SCF) iterative method developed by Hachisu (1986). In 
this method, the radial profile of the angular velocity or, equivalently, the specific angular 
momentum j — r^Q, is specified. Following Pickett et al. (1996), we fix j{r): 
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3 = 



JrO 



( 



Mr + Mc 



) 



(16) 



where jVo a-nd are constant, and is the disk mass within radius r. Setting q — 2 gives a 
j profile close to that used by Pickett et al. (1996). We begin the iteration with an arbitrary 
mass density p, from which we can calculate $s. From p and the above expression for j we 
also calculate Q (note that it still depends on the constant jro)- The relation (15) is then 
integrated to give the value of the enthalpy h: 



where the constant C and jro are determined from the boundary conditions p = at 
(r = Rin,z = 0) and at (r = Rout,z = 0). Here and Rout are the radial boundaries of 
the disk. The new density field is then calculated from h using the normalizing condition 
that Pmax = 1; which determines the polytropic constant K. Upon iterating this procedure, 
we converge to a model very close to equilibrium. 

The resulting disk model (with Rin = 0.25 and Rout = 1) has an Q profile close to 
Keplerian and a density profile displayed in figure 1. Note that the disk is rather thin, with 
an aspect ratio H/r varying between 0.1 and 0.2. 

As noted above, the ratio Mc/M^ is chosen in such a way that the Toomre Q parameter 
is initially close to unity. The radial profile of Q in the initial disk model is shown in figure 2. 
The minimum value of Q is approximately 1.1, and Q is close to unity over a large range of 
radii. We therefore expect strong non-axisymmetric gravitational instabilities to develop in 
this disk. 



Table 1 lists the parameters of the different runs we present below. Column 1 gives 
the label of the model. HD refers to a hydrodynamical run. Models T and P start with a 
purely toroidal and poloidal magnetic field, respectively. Column 2 gives the computational 
azimuthal domain and column 3 gives the highest fourier component of the gravitational 
potential. When rumax — 0, i.e. when only the m = component in the Fourier expansion 
of $5 is included, gravitational instabilities cannot develop (recall that Q > 1, so that the 
disk is stable against axisymmetric perturbations). Therefore, models Tl and PI enable us 
to study the evolution of MHD turbulence and to compare it with previous work and with 
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4. Results 
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Fig. 1. — Density contours in the (r, ;2)-plane of the initial disk model used in all the 
simulations. The contours shown are p = 10"^, 10-^ 10"^ 0.01, 0.1, 0.5, 0.7 and 0.9. The 
disk radial boundaries are Rin — 0.25 and Rout — 1- The central mass is twice that of the 
disk. 
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Fig. 2. — Radial profile of the Toomre Q parameter in the initial disk model. Q is close to 
unity over a large range of radii. 
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the 2D simulations of paper I. In models T2, T2iow, T2*, T3 and P2, both gravitational and 
magnetic instabilities develop. In model T2*, the non-axisymmetric part of $5 is included 
only after 6 orbits, i.e. after MHD turbulence has established itself. Column 4 gives the ratio 
of the volume-averaged thermal and magnetic pressures and column 5 gives the resolution 
{Nr, N^, N^) of the run. 

In all the models, an adiabatic equation of state is used. The computational domain 
extends radially from 0.1 to 1.4, and vertically from —0.2 to 0.2. In the azimuthal direc- 
tion, the computational domain extends from to either tt/2 or n. The smaller range is 
used in the rrimax = gravitationally stable cases. Indeed, Hawlcy (2000) and Papaloizou 
& Nelson (2003) have shown that an azimuthal domain of 7r/3 is generally sufficient to de- 
scribe the transport properties of MHD turbulence. When we allow for the development 
of gravitational instabilities, we restrict the azimuthal domain to the half disk [0, tt]. This 
saves computational time, but of course allows only even modes to develop. The focus of 
the paper is not on the detailed spectrum of modes which appear in a given disk model, 
however, but on the interaction between MHD turbulence and the largest scale gravitational 
modes. This interaction should not be particularly sensitive to whether an integer number 
of modes exactly fits in the half disk. 

Time is measured in units of the orbital period at the initial outer edge Rout = 1 of 
the disk model. Typical simulations are carried out for 8 to 10 orbits at this position. This 
corresponds to 60-80 orbits at the initial disk inner edge. The simulations are seeded by 
adding to the mass density at r > 0.4 random perturbations with a relative amplitude of 
5 X 10-3. 

We now describe in turn the hydrodynamical run, the simulations with only MHD 
turbulence, and the runs with both gravitational and magnetic instabilities. 

4.1. Control Hydrodynamical Run: Model HD 

The time evolution of the Fourier components of the density in the equatorial plane is 
shown in figure 3 for the modes m = 2, 4 and 6 (from top to bottom). The m — 2 mode grows 
at the beginning of the simulation and saturates after 4 orbits. Higher m modes emerge after 
about 3 orbits. Apart from the m = 4 mode, which may also be linearly unstable, the m > 2 
modes appear to be non-linearly excited. 

The development of a m = 2 spiral structure may be seen in figure 4, which shows the 
logarithm of the density in the equatorial plane at t = 4.27. (Note that the result of the 
simulation has been extended by symmetry to cover the range [0, 27r]). This mode is clearly 
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^ For this run, rrimax = when t G [0, 5.8], while rrimax = 8 for f > 5.8. 

Table 1: Model parameters. Column 2 gives the computational azimuthal domain, column 3 
gives the highest Fourier component of the gravitational potential included in the calculation, 
column 4 gives the ratio of the volume averaged thermal and magnetic pressures and column 5 
gives the resolution {Nr, N^, N^) of the run. Model HD is hydrodynamical. Models T and P 
start with a purely toroidal and poloidal magnetic field, respectively. When mmax — 0, the 
disk self-gravitating potential is forced to stay axisymmetric. In model T2*, rrimax = at 
the beginning of the run and is set to 8 after a few orbits. 
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Fig. 3. — Time evolution of the Fourier components of the density in the equatorial plane for 
model HD. The y-axis represents the ratio of the amplitude of the m-th Fourier component 
of the perturbed density to the unperturbed density. From top to bottom, the different curves 
correspond to the modes m = 2 {solid line), m = 4 {dashed line) and m = 6 {dotted- dashed 
line) . 
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global. Its pattern speed is Qp — 6.28, which means that corotation (the radius where the 
gas angular velocity matches the pattern speed) is located at the initial outer edge of the 
disk. Such a mode is predicted to emerge by linear stability analyses of self-gravitating disks 
(Papaloizou & Savonije 1991). The instability is due to the interaction between waves that 
propagate near the outer boundary and waves that reside inside the inner Lindblad resonance 
(where the pattern speed in the frame corotating with the planet matches the gas epicyclic 
frequency). This is located at r ~ 0.6 in our disk model. 

During the simulation, matter is driven toward the disk center by the gravitational 
torque associated with the spiral arms and, at t ~ 8, Q has become sufficiently high (^2) 
that the disk settles into a stable state. The results of this simulation are in agreement with 
theoretical expectations and with previous work, and show that the Poisson solver performs 
satisfactorily in the hydrodynamical regime. 

4.2. MHD Simulations in an Axisymmetric Gravitational Potential 

In the presence of a weak magnetic field, we expect our disk model to be unstable to the 
MRI, regardless of the field geometry. We first perform simulations in which only the MRI 
develops (models Tl and PI). This allows the properties of the ensuing MHD turbulence to 
be quantified and compared with previous work. To prevent the growth of non-axisymmetric 
gravitational instabilities, we retain only the m — component in the Fourier expansion of 

The resolution is (A^^, Nr^, Nz) — (128, 32, 128) and the azimuthal domain extends from 
to 7r/2, which would be equivalent to a resolution of 128^ over a range of 27r. 

4-2.1. Initial toroidal field: Model Tl 

We add to the equihbrium disk model described above a, {P) — S toroidal magnetic field 
and run the simulation for about 10 orbits (about 80 orbits at the initial disk inner edge). 

Figure 5 shows the time evolution of the volume-averaged Maxwell and Reynolds stress 
tensors and the corresponding a parameter (see eq. [14]). The Maxwell stress increases 
during the hnear phase of the instability. It then saturates after 4 orbits, when the MRI 
breaks down into turbulence. The presence of turbulence is seen in figure 6, which shows 
the density perturbation in the equatorial plane at t = 6.4. Turbulent fiuctuations are 
present over the full extent of the disk. It is clear from figure 5 that the Reynolds stress is 
significantly smaller than the Maxwell stress over the course of the simulation. This is in 
agreement with previous non self-gravitating global simulations of the MRI (Hawley 2000, 
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2001; Steinacker & Papaloizou 2002). The right panel of figure 5 shows the radial profile of 
a at the end of the simulation, i.e. at t ~ 10. The typical value of a is a few times 10~^, 
similar to what was found in previous simulations starting with a toroidal field with a net 
fiux (Steinacker & Papaloizou 2002). 

The Maxwell stress stays roughly constant during our simulation. This indicates that the 
resolution (128, 32, 128) is large enough for the turbulence to be sustained over the duration 
of the run. We therefore adopt it in the following runs (which of necessity are limited in 
time by the fact that mass is accreted onto the central mass). 



To investigate the sensitivity of our results to the initial field geometry, we run the same 
calculation as in model Tl but with an initial poloidal magnetic field. We calculate the field 
from the (toroidal) component of the vector potential in the initial disk model: 



This corresponds to 4 magnetic loops confined inside the disk. The first 2D simulations of 
a disk permeated by a weak vertical field (Hawley & Balbus 1991) showed the development 
and growth of "channel" solutions. In 3D, these solutions still exist but they quickly break 
down into turbulence, as predicted by the analysis of Goodman & Xu (1994). Turbulence 
is more rapidly established when the field varies on a fairly small scale, which motivates the 
above choice of A^. 

The radial and vertical components of the magnetic field are computed from and 
normalized such as to obtain the desired initial value of Since the linear growth of the 
vertical field is much more rapid than that of the toroidal field (see below) , we chose a much 
larger initial value of = 300. 

The properties of the turbulence arc similar to those found when the initial field is 
toroidal. Figure 7 shows the time evolution of the Maxwell stress for both models Tl and 
PI. As expected, the linear instability is much more vigorous when a vertical field is present, 
because of the growth of the channel solutions. However, in both cases the stress saturates 
at a similar value and the level of turbulence is comparable. The evolution of the Maxwell 
stress in PI is somewhat similar to what was obtained in paper I. The important difference 
is that in 3D, the stress saturates when turbulence is estabhshed and does not decay with 
time, as it does in 2D. 



Jf..2.2. Initial poloidal field: Model PI 




(18) 
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Fig. 4. — Logarithm of the density in the equatorial plane for model HD at t = 4.27. The 
linearly unstable two-arm global mode has become non-linear. It drives angular momentum 
outward and matter toward the central point mass. 




Fig. 5. — Left panel: Time evolution of the volume averaged Maxwell {solid line) and 
Reynolds {dashed line) stress tensors for model Tl. The Maxwell stress increases during 
the linear growth of the MRI (first 4 orbits). It then saturates when the instability breaks 
down into turbulence and stays roughly constant. At all time, the Reynolds stress is much 
smaller than its magnetic counterpart. Right panel: a vs. r at the end of the simulation, i.e. 
at t ~ 10. The typical value of a is a few times 10~^. 
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Fig. 6. — Density perturbation in the equatorial plane for run Tl at t = 6.4. Turbulent 
fluctuations are present over the whole extent of the disk. 
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Fig. 7. — Time history of the volume averaged Maxwell stress tensor for runs PI (dashed 
line) and Tl (dotted line). The linear instability is more vigorous when a vertical field is 
present, but the level of turbulence is similar in both cases. 
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4.3. Full MHD simulations 

In this section we report the resuhs of full 3D simulations including the develop- 
ment of both self-gravitational and magnetic instabilities. The resolution is {Nr, N^p, N^) = 
(128, 64, 128) and the azimuthal domain extends from to tt. 

4.3.1. Initial Toroidal Field: Models T2, T2* and T3 

In this sequence of models, we observe the simultaneous appearance of both MHD 
turbulence and the m = 2 spiral arm familiar from the hydrodynamical calculation. To 
better understand how angular momentum is transported in the disk, we compare the time 
evolution of the different stresses with those obtained in the models described in the previous 
section. 

Figure 8 shows the gravitational stress (T^™^) as a function of time for both models T2 
and HD. Somewhat surprisingly, the presence of both gravitational and MHD instabilities 
leads to an average (T^™^) reduced by a factor ~ 2 compared with the values obtained 
without a magnetic field. The magnetic torques do not lead to more vigorous gravitational 
instability. One possible explanation may be that turbulent motions tend to broaden the 
spiral arms by adding an extra fluctuating component to the thermal pressure, but there 
is more going on just this. Figure 8 also shows that the gravitational stress varies nearly 
periodically with time and can reach very small values. This behavior is in fact associated 
with the near disappearance of the spiral arms, as can be seen in figure 9. These snapshots 
correspond to a maximum and a minimum of the gravitational stress, respectively. The 
spiral arms are sharp at t = 4.95, whereas they lack definition at t = 5.09. The arms form, 
disperse, and reform. This periodic variation is also seen in figure 10, which shows the mass 
accretion rate onto the central mass as a function of time. As expected, the accretion rate 
has a periodic component with the same frequency as that found in the gravitational stress. 
The period in both cases is ~ 0.28. 

In figure 11, we compare the evolution of the Maxwell stress in models T2 and T2* 
(for which only the first 6 orbits, during which nimax — 0, are plotted). In contrast to the 
gravitational stress, the Maxwell stress is significantly larger when the disk is gravitationally 
unstable. This appears to be due to the systematic compression of the magnetic field lines 
along the spiral arms, as opposed, say, to an increase of the level of turbulent fluctuations. 
When the gravitational instability disappears after about 7 orbits in model T2, for example, 
the Maxwell stress decreases to the same value as in model T2*. 

As in the hydro model HD, the Toomre Q parameter rises throughout the body of the 
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Fig. 8. — Time evolution of the volume averaged gravitational stress tensor {T^^^'") for models 
HD {solid line) and T2 {dotted line). The level of transport by gravitational instabilities is 
significantly reduced when MHD turbulence is present. In model T2, the gravitational stress 
also varies periodically. 
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Fig. 9. — Logarithm of the density in the equatorial plane for run T2 at t = 4.95 {left panel) 
and t = 5.09 {right panel). These snapshots correspond to a maximum and a minimum of the 
gravitational stress, respectively. The spiral arms are sharp and clear at t = 4.95, whereas 
they appear blurred at t = 5.09. 




Fig. 10. — Time evolution of the mass accretion rate (mass per unit of time) onto the central 
mass for run T2. As expected, the accretion rate oscillates with the same period as the 
gravitational stress. 
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disk over the course of the simulation (as mass is transported toward the inner region), 
until the gravitational instability ceases. But even by t ~ 8, when there is no longer any 
gravitational transport, Q is still larger in model HD than in model T2. This is because 
the gravitational instability is stronger in the hydro dynamical case, and the disk is depleted 
more rapidly. 

To check the sensitivity of these results to our choice of initial conditions, we conducted 
the following experiment. In model T2*, the input parameters are the same as in model T2, 
but the non-axisymmetric part of $s is included only after 6 orbits, i.e. only after MHD 
turbulence has been firmly established. Figure 12 shows the time evolution of the volume 
averaged Maxwell and gravitational stress tensors for run T2*. Until t = 7, the development 
of MHD turbulence is the same as in model Tl. However, in the time interval t — 7-8, i.e. 
after gravitational instabilities have developed, {T^"-^) decreases to reach about one third 
of its value aXt — 1. The reason of this dechne is not completely clear. One possibility may 
be that the compression of the (randomized) magnetic field in the spiral arms leads to more 
efficient reconncction of the field lines. Another possibility is that gravitational stresses feed 
off the density fiuctuations generated by the MRl, thereby indirectly coupling the magnetic 
and gravitational energies. In any case, this behavior stands in contrast with was observed in 
run T2, where gravitational instabilities developed while the magnetic field was still ordered. 
In model T2*, the gravitational stress tensor is roughly a factor of 2 smaller than in model 
T2, but shows the same periodic variations. 

Our next comparison run, model T3, differs from model T2 only in the number nimax 
of fourier coefficients in the expansion of $s- (T3 has rrimax = 16, T2 has rrimax = 8.) Once 
again, very similar results emerge, and the choice of rrimax does not appear to be critical (cf. 
§ 4.3.3). 

To summarize: the evolution of a purely toroidal field in a gravitationally unstable disk 
leads to a reduction and strong periodicity in the gravitational stress (compared with a purely 
hydro dynamical model). Compared with gravitationally stable models, the Maxwell stress is 
larger or smaller depending repsectively on whether gravitational instabilities develop at the 
same time as MHD turbulence (magnetic field alignment in spiral arms) or after turbulence 
is established (reduction in magnetic stress as gravitational stress develops). The behavior 
of an initial poloidal field is considered next. 
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Fig. 11. — Time evolution of the volume averaged Maxwell stress tensor for runs T2* [dashed 
line), before gravitational transport is turned on, and T2 (dotted line). The Maxwell stress 
is larger when the disk is gravitationally unstable. This appears to be due to a compression 
of the magnetic field lines along the spiral arms. When the disk in model T2 becomes 
gravitationally stable (after i — 7), the Maxwell stress decreases down to the same value as 
in model T2*. 
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Fig. 12. — Time evolution of the volume averaged Maxwell {solid line) and gravitational 
{dashed line) stress tensors for run T2*. In this run, the non-axisymmetric part of $s is 
included only after 6 orbits, i.e. only after MHD turbulence has been firmly established. 
The development of gravitational instabilities coincides with a significant decrease of the 
Maxwell stress tensor. This may be due to reconnection of the (randomly oriented) field 
lines in the spiral arms. The gravitational stress tensor has the same amplitude and periodic 
variations as in run T2. 
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4-3.2. Initial Poloidal Field: Model P2 

Do our toroidal field findings extend to poloidal field behavior? To answer this question, 
we begin with an initial poloidal field, calculated as in section 4.2.2 above. Again, we start 
with {(3) = 300. Except for the initial field geometry, model P2 is the same as model T2. 

Figure 13 shows the evolution of the gravitational stress tensor for both models P2 and 
HD (this is the equivalent of figure 8). As in the case of a toroidal field, a non-axisymmetric 
m = 2 spiral grows and becomes nonlinear in model P2. The fact that gravitational insta- 
bilities develop earher in model P2 than in models T2 and HD (see figures 8 and 13) appears 
to be due to the fact that the strong hnear magnetic instabihty associated with the poloidal 
field produces large perturbations of the density. Once again, we find that the gravitational 
stress is smaller than in model HD, and varies periodically with time with a period ~ 0.38 
somewhat larger than for model P2. 

Figure 14 shows the evolution of the Maxwell stress for models PI and P2 (this is 
the equivalent of figure 11). Once again, the Maxwell stress is larger when the disk is 
gravitationally unstable and gravitational instabilities develop at the beginning of the run 
(model P2). 

Since the state of MHD turbulence in a saturated disk is the same whether an initial 
poloidal or toroidal field is used, we have not run a case "P2*" with the non-axisymmetric 
part of $5 added later. Such a run is expected to be very similar to T2*, since the initial 
turbulent states are similar. 

The history of run P2 and the above argument together suggest that toroidal and 
poloidal initial fields behave very similarly. 

4-3.3. Modal analysis 

The full MHD simulations described above suggest that a general feature of the evolution 
of gravitationally unstable turbulent disks is a periodic modulation of the gravitational stress. 
To understand the reason for this modulation, we now analyse in more detail the unstable 
modes that appear in models HD, T2 and P2. We are in particular interested in determining 

the power spectrum as a function of mode frequency a. Following Papaloizou & Savonije 
(1991), we Fourier transform in time at each radial zone r and at a fixed azimuth 0o the 
function p{r,(j)Q,t). To get the spectral time evolution, we carry out each Fourier transform 
over a series of 4 distinct time intervals. The number of time-steps used in each time interval 
gives a finite frequency resolution da/27i = 0.3. The contours of constant power are then 



-22- 



1.5x10-^ 



o 

c 1.0x10-^ 

CD 



5.0x10- 



-5.0x10-5 








4 
t/To, 



Fig. 13. — Time evolution of the volume averaged gravitational stress tensor for models HD 
{solid line) and P2 (dotted line). The gravitational stress in model P2 is smaller than in 
model HD and varies periodically with time. 
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Fig. 14. — Time evolution of the volume averaged Maxwell stress tensor for models PI 
{dashed line) and P2 {dotted line). The Maxwell stress is larger when the disk is gravitation- 
ally unstable and gravitational instabilities develop at the beginning of the run. 
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plotted as a function of frequency and radius for the various time intervals. 

Figure 15 shows the contours for model HD. The different panels correspond to different 
time intervals. In the first panel, i.e. at time t ~ 2.4, we see the presence of a mode with 
frequency (t/2t: — 2.5 which extends over the whole disk. Its amplitude is small, peaking 
at about 3 x 10~^. This mode is still present in the second panel, at time t ~ 3.4, with 
a similar amplitude structure. However, a second low frequency mode with a /2tt = 1 (i.e. 
with a corotation radius at the disk initial outer edge) is now apparent. Its amplitude is 
significant only in the disk outer parts, where it peaks at 0.2. This mode subsequently grows 
and completely dominates the high frequency {a/2'ri = 2.5) mode in the third panel, i.e. at 
t ~ 4.4. There its amplitude peaks at 0.7. This mode can of course be identified with the 
two-arm spiral seen in figure 4. At t — 4.4 the mode is non-linear. Its amplitude does not 
increase further, as can be seen in the fourth and last panel, i.e. at t ~ 5.4. At this later 
time, only the frequency of the mode has changed. It is now around 1.5. Both the finite 
frequency resolution da and the increase of the central mass due to accretion may account for 
this shift in frequency. A third mode with a frequency twice that of the low frequency mode 
is seen in the last two panels. The relationship between the frequencies of these two modes 
and the fact that their radial structure is very similar suggests that they are harmonics of 
each other. 

Figure 16 shows the contour plots for model T2. In the first panel, at t ~ 2.4, there 
is no dominant mode. Instead, there is a large number of high frequency (with mostly 
(T/27r = 1.3-3.8) perturbations. The amplitude of these fiuctuations is a few times 10~^. 
They are associated with the growing MRI. In the second panel, at t ~ 3.4, two modes 
emerge, but their amplitude is still rather low. However, these modes subsequently grow and 
are clearly seen with a larger amplitude in the third panel, at i ~ 4.4. One of these modes has 
a frequency a/2Tr — 1 and is the same as that seen in the hydrodynamical simulations. Its 
corotation radius is located at the disk initial outer edge. Its amplitude peaks at a value of 
about 0.5 in the outer parts of the disk. The other mode has a frequency a/27r = 2.5 and an 
amplitude ~ 0.2 constant over the whole disk. In particular, its amplitude in the disk inner 
parts is larger than that of the other mode. This mode is probably the same as that seen 
with a lower amplitude at early times in the hydrodynamical simulations (panels 1 and 2 of 
figure 15). This suggests that this mode is a disk eigenmode which is excited in model T2 by 
the high frequency motions associated with the turbulence. We expect nonlinear coupling 
between these two modes to give rise to beat oscillations, i.e. oscillations with a frequency 
being a hnear combination of the frequencies of the two modes. We have noted above that 
the gravitational stress tensor in model T2 oscillates with a period ~ 0.28. This is consistent 
with the frequency of the oscillations being a/2n ~ m((THF — ctbf) = 3, where cthf and (Tbf 
are the frequencies of the high and low-frequency modes, respectively. This suggests that 
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Fig. 15. — Contours of constant Fourier power as a function of dimensionless mode frequency 
o jl-n (vertical axis) and radius r {horizontal axis) for model HD. The different panels corre- 
spond to the time intervals 1.94-2.91 {upper left), 2.91-3.88 {upper right), 3.88-4.85 {lower 
left) and 4.85-5.83 {lower right). The contour levels in each of these panels are, in unit 10~^: 
(0.22; 0.29; 0.44; 0.87; 1.1; 1.6; 3.2), (2.1; 2.6; 3.4; 5.2; 10; 15), (8.4; 11; 14; 21; 42) and (9.4; 
12; 16; 23; 47; 62), respectively. One mode (a two-arm spiral) dominates the spectrum in 
the non-linear stage (third panel), with a frequency (t/2t: — 1. 
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the oscillations of the gravitational stress tensor result from a nonlinear coupling between 
these two modes. 

Figure 17 shows the contour plots for model P2 in the time interval 2.72-3.88. the 
situation is similar to model T2. Here again two modes of comparable amplitude are present. 
One of this mode has a frequency (7/27r = 1 and can be identified with the mode which 
emerges in the hydrodynamical simulations. The other mode has a frequency (7/27r = 2.1. 
Given the finite frequency resolution dojln^ this second mode may be the same as that 
identified in model T2. However, it may also be a different mode with a lower frequency. In 
any case, the situation is qualitatively the same here as in model T2. Since the two modes 
have comparable amplitude, they interact nonlinearly, which results in a periodic modulation 
of the gravitational stress. 

We now vary some of the parameters of model T2 in order to examine the sensitivity of 
these physical results to the numerical input parameters. Figure 18 shows the contour plots 
for model T2;o«, in the time interval 2.72-3.88. The resolution for this run in the radial and 
vertical directions is half that of model T2. The similarity between this plot and the third 
panel of figure 16 demonstrates that our results do not depend strongly on the numerical 
resolution. Figure 19 is the same as figure 18 but for model T3. Again, it is very similar to 
the third panel of figure 16. Since, in model T3, only the parameter Vfimax is different from 
model T2, figure 19 suggests that the limited number of Fourier components inchided in the 
calculation of the self-gravitating potential does not qualitatively affect the main physical 
results presented in this paper. 

In conclusion, models P2, T2/oui and T3, taken together, suggest that the physical results 
described in this paper are insensitive to both the numerical setup of the simulations and 
the initial magnetic field topology. 

5. Discussion 

In this paper, we have presented the first 3D numerical simulations of the evolution of 
self-gravitating and magnetized disks. We have investigated disks in which only gravitational 
or magnetic instabilities develop, and disks in which both types of instabilities occur. 

When no magnetic field is present, self-gravitating disks are unstable when the Toomre 
Q parameter is close to unity. The spectrum of unstable modes in that case is dominated 
by a large scale two-arm spiral whose corotation radius is located near the disk outer edge. 
The instability is due to the interaction between waves which propagate near the disk outer 
boundary and inside the inner Lindblad resonance (ILR), respectively (Papaloizou & Savonije 
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Fig. 16. — Same as figure 15 but for model T2. The contour levels in each panel are the 
same as in figure 15. Here two modes dominate the spectrum in the nonlinear phase (third 
panel), with frequencies a/2T^ — 1 and (7/27r = 2.5, respectively. 
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Fig. 17. — Same as figure 15 but for model P2 and only in the time interval 2.72-3.88. The 
contour levels used here are, in unit 10~^: (5.1; 6.4; 8.5; 13; 25). 
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Fig. 18. — Same as figure 15 but for model T2;ou, and only in the time interval 2.72-3.88. 
The contour levels used here are, in unit 10~^: (5.1; 6.4; 8.5; 13; 25). 
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1991). 

When a magnetic field is present, two large scale modes grow in the disk. They both 
have m — 2. One of this mode is the same as that seen in the hydrodynamical simulations. 
The second mode has a higher frequency. The ILR of the low-frequency mode, which is at 
r ~ 0.6, is very close to the corotation radius of the high-frequency mode, which is at r ~ 0.5. 
Such a pair of modes was seen in the hydrodynamical calculations of self-gravitating disks 
performed by Papaloizou & Savonije (1991). There, both modes were unstable because of 
the particular vortensity profile. In our simulations, in the absence of magnetic field, only the 
low-frcqucncy mode is unstable. The high-frequency mode seems to be part of the spectrum 
of normal modes in that case, but it does not grow. The presence of MHD turbulence does 
not modify the spectrum of large scale (comparable to the disk radius) modes, as it acts 
on scales hmited by the disk thickness. However, the fact that the high-frequency mode is 
unstable in the MHD simulations suggests that turbulence acts as a source for high-frequency 
oscillations. 

Nonlinear coupling between these two modes leads to an oscillation of the gravitational 
stress tensor. Note that such a coupling between two modes with coinciding resonances has 
been suggested to explain some of the features seen in numerical simulations of galactic disks 
by Tagger et al. (1987). These authors argued that the proximity of the resonances made 
the coupling very efficient. The oscillation of the gravitational stress tensor is accompagnied 
by the periodic disappearance of the spiral arms in the disk. Also, the peak value of this 
stress is decreased by about half compared to the hydrodynamical simulations. 

The results reported here are robust and do not depend on the geometry of the magnetic 
field. They have important consequences for disks around AGN and protoplanetary disks. 
They first show that accretion of a self-gravitating disk onto the central star is slowed down 
when a magnetic field is present. They also show that the accretion is time-dependent, with 
a characteristic timescale for the variability being on the order of a fraction of the dynamical 
timescale at the outer edge of the region where the instabihties develop. 

As mentioned in the introduction, protoplanetary disks arc probably self-gravitating in 
the early phases of their evolution. For a disk of about 100 AU, the work presented here 
suggests variabihty on a timescale ~ 10^ years. The periodicity in the spatial distribution of 
knots in jets emanating from such objects is in the range 10-10^ years (Reipurth 2000), and 
is usually thought of as being produced by a time-dependent accretion in the central parts 
of the disk. The simulations presented in this paper suggest that periodic modulations of 
the accretion rate might well be the result of the interplay between gravitational instabilities 
and MHD turbulence, a far from obvious source. Note that the first detection of near-lR 
variability in a sample of Class I protostars was performed recently by Park & Kenyon (2002). 
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However, the poor time coverage of their data prevents a useful measure of the variabihty 
timescales to be extracted. 

Disks around AGN display time-dependent phenomena on a large range of timescales 
(Ulrich et al. 1997). The dynamical timescale for a disk orbiting a 10^ solar masses black 
hole at 10~^ parsecs is 9.3 years, and variations are observed on timescales up to years. This 
again is consistent with the processes described in this paper. 
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Fig. 19. — Same as figure 15 but for model T3 and only in the time interval 2.72-3.88. The 
contour levels used here are, in unit 10"^: (5.1; 6.4; 8.5; 13; 25). 
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